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1 Summary 

Small biomolecular circuits with two distinct and stable steady states have been identified as essential 
components in a wide range of biological networks, with a variety of mechanisms and topologies giving 
rise to their important bistable property. Understanding the differences between circuit implementations 
is an important question, particularly for the synthetic biologist faced with the challenge of determining 
which bistable circuit design out of many is best for their specific application. In this work we explore 
the applicability of Sturm's theorem — a tool from 19th-century real algebraic geometry — to comparing 
"functionally equivalent" bistable circuits without the need for numerical simulation. We consider two 
genetic toggle variants and two different positive feedback circuits, and show how specific topological 
properties present in each type of circuit can serve to increase the size of their operational range. The 
demonstrated predictive power and ease of use of Sturm's theorem suggests that algebraic geometric 
techniques may be underutilized in biomolecular circuit analysis. 

2 Key words 

synthetic biology; biological circuit; bistability; algebraic geometry; Sturm's theorem 

3 Introduction 

The field of synthetic biology has rapidly matured to the point where it is now possible to produce complex 
synthetic networks with prescribed functions and level of performance [1]. As in other fields of engineering, 
advances have been enabled by the use of small interchangeable modules that are "functionally equivalent" 
from an input-output perspective [2]. Bistable circuits — which play a role in essential biological processes 
including cell fate specification [3], cell cycle progression [4], and apoptosis [5] — make up a particularly 
large and diverse functionally equivalent set [6]. Effectively characterizing and comparing these biocircuits 
is crucial for determining which architecture is in some sense optimal for a particular context. 

Ordinary differential equation (ODE) models can be powerful tools for contrasting different biocircuits' 
"dynamic phenotypes" (see, e.g., [7]); however, as circuit size increases, the usefulness of such models can be 
limited by their complexity. Many of the relevant parameters are often unknown, and while computational 
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techniques have advantages, analytical criteria that focus on topology can provide a more exact assessment 
of a module's properties [8,9]. A novel analytical tool that can provide topology-based insights can be 
found in Sturm's theorem [10], developed in 1835 as a solution to the problem of finding the number of 
real roots of an algebraic equation with real coefficients over a given interval. Despite its predictive power, 
this "gem of 19th century algebra and one of the greatest discoveries in the theory of polynomials" [11] 
remains unexploited as a tool for synthetic biology. 

In this work we demonstrate an approach to bistable circuit discrimination based on Sturm's theorem 
that gives boundaries of the regions of bistability as exact analytic expressions, eliminating the need 
for numerical simulation. We compare the regions of bistability for two variants of the classic double- 
negative toggle switch as well as two positive feedback circuits, one of which is based on the bacteriophage 
A promoter Prm- Overall our results demonstrate the usefulness of algebraic geometric techniques for 
studying functionally equivalent bistable circuits and suggest a new and simple-to-implement method for 
choosing between them for synthetic biological applications. 

4 Mathematical preliminaries 

4.1 Sturm's theorem 

Sturm's theorem gives the number of distinct real roots of a univariate polynomial f(x) in a particular 
interval. To apply the theorem, we must first construct the Sturm sequence, a set of polynomials J- = 
{/o,/i, • • -,fm} defined as: 

/o = /(*), 

h = fix), 

h = -rem(/ 0 ,/i), 

h = -rem(/i,/ 2 ), 

f m = -rem(/ m _ 2 ,/m-i), 
0 = -rem(/ m _i,/ m ), 

where rem(/j, fi+i) is the remainder of the polynomial long division of fa by fi+i- The sequence ends at 
f m , when f m -i divided by f m gives a remainder of zero. For a polynomial of degree n, there are m < n + 1 
Sturm polynomials in the sequence. 

Theorem 1 (Sturm's theorem) Letf(x) be a real-valued univariate polynomial and a, b € RU{— oo, +00} , 
with a <b and f(a),f(b) ^ 0. Then the number of zeroes of f{x) in the interval (a,b) is the difference 

var{F, a) — var(F, b) , 

where J- is the Sturm sequence of f(x), and the variations war(J 7 , a) and var{J-,b) are the number of 
times that consecutive nonzero elements of the Sturm sequence — evaluated at a and b, respectively — have 
opposite signs. (Adapted from [12].) 

4.2 Number of steady states and bistability 

Our approach involves identifying regions of bistability by finding conditions that lead to three steady 
states, without requiring numerical determination of the exact values or stability of the equilibrium points. 
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While it is in general not possible to draw conclusions on the stability properties of equilibria by simply 
counting their number, the circuits under consideration enjoy two important properties — namely, they 
are dissipative and their linearizations are positive — that allow us to relate their degree and number of 
equilibria to the stability properties of each equilibrium (see Supporting Information). For such circuits, 
when three equilibria are present, two of them must be stable and one must be unstable. 



5 Results 



5.1 Genetic toggle circuit comparison 

A recent study identified a set of eleven minimal bistable networks (MBNs), simple two-gene circuits with 
the capacity for bistability that do not also contain a smaller bistable subcircuit [13]. One of these MBNs, 
consisting of two mutually-repressing dimeric proteins (Fig. 1A, top), resembles the toggle that was among 
the very first synthetic biocircuits [14]. Another is a toggle variant in which one of the repressors functions 
as a monomer (Fig. 1A, bottom). To our knowledge no monomeric repressor yet exists; however, given the 
protein engineering capabilities of synthetic biology, predicting how such a toggle will perform relative to 
a more natural functional equivalent is worthwhile. We refer to these two MBNs as the dimer-dimer (DD) 
and monomer-dimer (MD) toggles, respectively. 
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Figure 1: (A) Dimer-dimer (top) and monomer-dimer (bottom) toggle switches. (B) Bistable regions for 
the monomer-dimer and dimer-dimer toggles. 

Beginning with a chemical reaction network formulation and assuming mass-action kinetics we derive ODE 
models of the two toggles (Eqs. (SI) and (S2)). At equilibrium the concentrations of Pi and P 2 in the 
MD system are given by 

P leq = (1 + (P 2eq /K 2 ) 2 y 1 l3 1 X ltot , P 2eq = (1 + {Pxe q lK md )Y X hX 2tot , (1) 

and in the DD system, 

P leq = (1 + {P 2eq /K 2 f) ~ l ^X Uot , P 2eq = (1 + {P leq /K dd ) 2 ) ' X hX 2tot , (2) 
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where X{ tot is the total amount of gene i, is the ratio of basal production rate to degradation rate for 
protein i, and K m ^ K^, and K2 are Michaelis constants (see Supporting Information for details). We 
may nondimensionalize these equations by scaling the protein concentrations relative to their respective 
Michaelis constants, and the DNA concentrations accordingly: X 2tot = P2(X2 tot /K 2 ), P2 eq = P2e q /K 2 , 
Pieq = Pieq/Kxd, and X ltot = j3i(X Uot / K xd ) . Systems (1) and (2) can then be written as polynomials in 

Pleq'- 

^3 -~-2 2 

fmd(Pleq) = Pleq ~ {X\tot ~ tyPleq ~ (%X\ tat ~ X 2tot ~ l)Pl eq ~ X ltot = 0 (3) 

and 

fdd{Pleq) — Pleq ~ Xl tot Pl eq + 2P leq - 2X ltot P leq + (X 2tot + l)Pl e q ~ X Uot = 0 , (4) 

for the MD and DD toggles, respectively. Every positive root of these equilibrium polynomials gives 
a positive steady state concentration for every other circuit component as well. To find the regions of 
bistability in the plane of X\ tot and X 2to ty we construct the Sturm sequence T associated with each circuit's 
equilibrium polynomial (Eq. (3) or (4)), evaluate T at P\ eq — > 0 and P\ eq — > +oo, and find the conditions 
leading to a variation difference var{F, 0) — var(J 7 , +oo) = 3. (Sturm sequences are given in the Supporting 
Information.) 

The MD toggle Sturm sequence J- m d nas a maximum possible variation of 3 and only one combination of 
inequalities that can give rise to bistability: when var(F m d-,Q) = 3 and var(Fmd, +oo) = 0. In contrast, 
the DD toggle sequence T^d could in principle yield five or four positive steady states; however, only three 
are admitted as there are no combinations of inequalities that have a variation difference of 5 or 4 and 
are logically consistent. All possible sign sets for the DD toggle are listed along with their allowabilities 
in Tables SI and S2. The analytic expressions for the two regions of bistability are Eq. (S14) (MD toggle) 
and the intersection of Eqs. (S15) and (S16) (DD toggle). We find that the DD toggle operates over a 
substantially greater range than the MD toggle (Fig. IB), indicating that the DD architecture is more 
functionally robust to variations in circuit DNA concentration (given fixed rate parameters). Furthermore, 
the DD switch can operate with significantly lower concentrations of DNA: a >50% reduction in X 2 tot an d 
>75% reduction in X Uot . 

Recognizing that some of the mathematical tools used may be unfamiliar, some computational validation 
of our results may be of value. Computational tests and results that perfectly match the Sturm results are 
described in the Supporting Information. 

5.2 Single gene circuit bistability 

The single gene system consisting of bacteriophage A repressor and its promoter Prm (with its three 
operator sites OR1, OR2, and OR3) also exhibits bistability [15]. We can compare the bistability region 
of this multi-operator circuit with that of a simple positive feedback circuit containing only one operator 
site. 

A dimensionless model of the A single gene system is given in [15]. At steady state the concentration of 
protein is given by 

Wi^Pjg + 7^e 9 - aaxP^ + 7 P e 3 g - P^ + - 1 = 0 , (5) 

where 7 is the rescaled degradation rate constant, a represents the increase in protein production resulting 
from dimer binding to OR2, and a± and 02 are the relative (to OR1) affinities for OR2 and the negatively- 
regulating OR3, respectively. (For simplicity we set the gene copy number equal to one.) With a\ = 2 and 
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o 2 = 0.08 [15], the associated Sturm sequence Jr> RM has only two sign sets with var^p^, 0) = 5 and one 
set with var (-7"p RM , +oo) = 2 that are logically consistent and together give bistability (Table S3). 

In contrast, for a single gene positive feedback system with one operator site for its dimeric protein (MBN 
kqw in [13]), rescaled as in Eq. (5), we have: 

7 P e 3 ? - aP 2 eq + 7 P eq - 1 = 0 . (6) 

As with the MD toggle, this polynomial also has a maximum possible variation of 3 and thus only a single 
combination of inequalities that give rise to bistability, in the region given by Eq. (S17). 

The bistable regions (in a-7 space) for these single gene systems are shown in Fig. 2A. It can be seen 
that the A repressor circuit functions over a larger range and with lower values of the degradation rate 
constant. Interestingly, with a = 11 ( [15] and references therein) the single operator circuit would only 
just barely function as a bistable circuit, and any small fluctuation in circuit parameters would render it 
nonfunctional. 




Figure 2: (A) Bistable regions for the single operator positive feedback circuit and the A repressor-PRM 
system with relative affinities a\ = 2 and 02 = 0.08. (B) Region of bistability for the A-Prm system as a 
function of protein production enhancement a, protein degradation rate 7, and relative affinity for OR3 

Sturm's theorem may also be used to determine how the strength of the negative feedback ((72) affects 
bistability. Keeping o\ = 2, and using a = 11 and 7 = 4.5 (centered in the bistable region at a = 11; see 
Fig. 2A), we find that 02 can increase twelve-fold to ~ 0.96 before bistability is lost. In general, significant 
increases in a 2 require similar increases in a for bistability to maintained, with the range of allowable 7 
narrowing as result (Fig. 2B). 

6 Discussion 

We have shown how real algebraic geometry, specifically Sturm's theorem, can be used to compare and con- 
trast functionally equivalent bistable biocircuits. We take advantage of the fact that the circuits discussed 
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here, like many in biology, are dissipative and positive in their linearization, which allows us to ascertain 
the stability of their equilibria without computation. For each circuit the bistability region is given as a 
set of exact inequalities, simplified by the combination of the (often unknown and fixed) rate parameters 
into rescaling factors. Though applied to only four different bistable circuits here, this same approach may 
be used for any system that at steady state can be described by a univariate polynomial. 

Comparing two different genetic toggle variants, we see that one consisting of two dimeric repressor species 
has a wider operational range than one composed of one dimeric and one monomeric repressor. This 
result demonstrates the general importance of comparing functionally equivalent circuits, with particular 
relevance for the design of larger circuits that rely on bistable modules: given the uncertainty that exists 
in parameter values and DNA concentration (when DNA is in the form of plasmids without strict copy- 
number control), a bistable circuit variant that works over a wide range of parameters would most likely 
be preferable. 

Our results also demonstrate the benefit of an additional operator site for a real single gene positive feedback 
system: if the A repressor-Pjjjvi system did not have both OR1 and OR2, the enhancement a = 11 would 
barely be sufficient for the system to be bistable. We also see that the negative feedback at OR3 is small 
enough to not significantly affect bistability. This suggests that the promoter architecture of the A system 
may have evolved to allow for both robust bistability due to the positive feedback as well as reduced 
variability or other benefit of negative autoregulation. 

We have used Sturm's theorem to distinguish between biocircuits known to exhibit bistability. However, 
it may also be used (like Chemical Reaction Network Theory, previously [13,16]) to predict new bistable 
topologies or rule out those circuits that do not have this capacity, since only those circuits with a variation 
difference var{J-, 0) — war (J 7 , +oo) = 3 for some sets of parameters can be bistable. Altogether the addition 
of real algebraic geometry to the general synthetic biology toolbox can only help improve the process of 
biological circuit design and analysis. 
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SI Chemical reaction networks 



With the exception of the model for the multiple-operator Prm promoter of bacteriophage A, the various 
circuits highlighted in the main text were initially predicted to exhibit bistability using a chemical reaction 
network (CRN)-based topological survey [1]. Each of these CRNs contains reactions representing basal 
protein production and degradation: 



Xi-^Xt+Pt, X 2 =*2* x 2 + P 2 , P x 



P, 



for genes Xj and proteins Pj. The other reactions that uniquely define each circuit are: 



k c p 

X2 + Pi s X2P1 

k c R 

2P 2 P 2 P 2 

kkR 
k n p 

X\ + P2P2 s X\P 2 P 2 

k n R 



Monomcr-dimer (MD) 
toggle 



2Pi Pi Pi 

Kr 

2P 2 ^ P 2 P 2 

kkR 

X\ + P 2 P 2 v XlP 2 P 2 

k n R 

X 2 + Pi Pi v X 2 PiPi 

k D R 



Dimer-dimer (DD) toggle 



Hie F 

2P 2 ^ P 2 P 2 
kkR 

kqF 

X 2 + P 2 P 2 s X 2 P 2 P 2 

kqR 

k w 

X 2 P 2 P 2 ^ X 2 P 2 P 2 + P2 



Single-operator 
positive feedback circuit 



PiPi represent dimeric species, and XiPj and XiPjPj represent monomers and dimers bound to the gene 
promoters. The various ODE sets were derived from these CRNs under the assumption of mass action 
kinetics and simplified using the fact that the total concentrations of each gene (in bound and unbound 
form) are conserved. 
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52 Ordinary differential equation models 

From the chemical reaction network formulation and assuming mass-action kinetics we can derive the 
following sets of ordinary differential equations (ODEs) that describe the circuit dynamics. For the MD 
toggle: 

Pi '(*) = -k degX P x (t) - k cF X 2 (t)-P x (t) + k basX X x (t) + k cR {X 2tot - X 2 (t)) (Sla) 
P 2 \t) = ~2k kF P 2 {t) 2 - k deg2 P 2 (t) + 2k kR P 2 P 2 (t) + k bas2 X 2 (t) (Sib) 
iW(t) = k kF P 2 (tf - k kR P 2 P 2 {t) - k nF P 2 P 2 (t)-X x (t) + k nR {X ltot - X x {t)) (Sic) 
X x '(t) = k nR {X Uot - X x {t)) - k nF P 2 P 2 {t)-X x {t) (Sid) 
X 2 '(t) = k cR (X 2tot ~ X 2 {t)) - k cF P 1 {t)-X 2 {t) , (Sle) 

the DD toggle: 

Pi '(*) = -2k lF P 1 (t) 2 - kdeglPl{t) + 2k iR P X P X (t) + kbaslXxit) (S2a) 

Pi Pi '(t) = fc iF Pi(t) 2 - k iR P x P x {t) - k oF P x P x (t)-X 2 (t) + k oR (X 2tot - X 2 (t)) (S2b) 

P 2 '(t) = -2k kF P 2 {tf - k deg2 P 2 {t) + 2k kR P 2 P 2 (t) + k bas2 X 2 {t) (S2c) 

P 2 P 2 'W = k kF P 2 (t) 2 - k kR P 2 P 2 {t) - k nF P 2 P 2 (t)-X x (t) + k nR {X Uot - X x (t)) (S2d) 

Xi '(t) = fc„ fl (X ltot - X x (t)) - fc„FP 2 P 2 (i)-Xi(<) (S2e) 

X 2 '(t) = k oR (X 2tot - X 2 {t)) - k oF P x P x (t)-X 2 (t) , (S2f) 

and single-operator positive feedback circuit: 

P 2 '(t) - fcba, 2 X 2 (t) - k deg2 P 2 (t) - 2k kF P 2 (t) 2 + 2k kR P 2 P 2 (t) + k w (X 2tot - X 2 (t)) (S3a) 

P 2 P 2 '(*) = fc fc FP 2 (t) 2 - k kR P 2 P 2 (t) - k qF P 2 P 2 (t)-X 2 {t) + k qR {X 2tot - X 2 (t)) (S3b) 

X 2 '(t)=A; gi i(X 2toi -X 2 (t))- k qF P 2 P 2 (t)-X 2 (t) , (S3c) 

where the variable names are as in [1]: Xj is the concentration of free (i.e., unbound by repressor) gene 
i, Xj tot is the total amount of X{ in the circuit (bound and unbound), and P{ and PiPi represent the 
monomeric and dimeric forms of protein i, respectively. The various k x are the reaction rates, and in the 
positive feedback circuit, k w > kb aS 2 is assumed. 

53 Derivation of the toggle circuit equilibrium polynomials 

Using (SI) and (S2), with the left-hand sides set equal to zero, we can derive the univariate equilibrium 
polynomials used in the application of Sturm's theorem. 

For the DD toggle, we subtract (S2d) from (S2e): 

0 = k nR {X Uot - X x ) - k nF P 2 P 2 -X 1 

- k kF Pl + k kR P 2 P 2 + k nF P 2 P 2 -X x - k nR (X ltot - Xi) 

=> P2P2 = ^P 2 2 = p- , (S4) 

where for simplicity of notation we use Xj, Pj, and PiPi to mean the equilibrium concentrations. We then 
plug this expression into (S2e) to get 

0 = k nR (X ltot - X x ) - k nF {P 2 2 /k kD )-X l 

= (X ltot - X x ) - (P 2 2 /(k kD k nD ))-X l 



Y _ k kD k nD Xi tot 

Xl ~ k kD k nD +P? ■ (S5) 
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Similarly we subtract (S2b) from (S2f) 



0 = k oR (X 2tot - X 2 ) - k oF P l P 1 -X 2 

- k lF Pf + k m P 1 P 1 + k oF P 1 P r X 2 - k oR (X 2tot - X 2 ) 

PlPl = ^Lpl = ( S6 ) 
kiR k iD 



and plug the resulting expression back into (S2f) to get 

0 = k oR (X 2tot - X 2 ) - k oF (P?/k iD )-X 2 
= (X 2tot - X 2 ) - (P?/(k iD k oD ))-X 2 

^ X2 ~ k lD k oD + P? ■ (S7) 

Substituting Eqs. (S4)-(S7) into (S2a) and (S2c) gives the equilibrium concentrations of P\ and P 2 in the 
DD toggle shown in the main text: 

P leq = (1 + (P^/K^y'PiXuot , P2 £q = (1 + {Px m lK dd fY X p 2 X 2tot , (S8) 

where = kj, as i / k^egi is the ratio of basal production rate to degradation rate for protein i, and the 
Michaelis constants K d d = (^d^od) 1 ^ 2 and K 2 = (kkDknD) 1 ^ 2 represent the protein concentrations that 
yield 50% of the maximum production rate of their respective targets. The combination of these two 
expressions gives the DD equilibrium polynomial. 



For the MD toggle, we have again that 



and 



but also (from Eq. (Sle)) that 



P2P2 = p- (S9) 
kkD 



kkDk n DX\ tot /om\ 



kkDKD + P 2 



Y k c pX 2tot 

X 2 = 7 —5- • (Oil) 

kcD + Pi 

Substituting Eqs. (S9)-(S11) into (Sla) and (Sib) gives the equilibrium concentrations of Pi and P 2 in the 
MD toggle shown in the main text: 

Pleq = (1 + (P2e q /K 2 ) 2 )- l ^X ltot , P 2eq = (l + (P Ug /K md )) ^ /3 2 X 2tat . (S12) 

As previously, is the ratio of basal production rate to degradation rate for protein i, and the Michaelis 
constants K rn d = k cF > and K 2 = {kkokno) 1 ^ 2 represent the protein concentrations that yield 50% of the 
maximum production rate of their respective targets. The combination of these two expressions gives the 
MD equilibrium polynomial. 
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S4 Sturm polynomials 

Sturm polynomials can be rather long and complicated functions; however, they can be easily generated 
using software capable of symbolic manipulation, e.g., Mathematica. We present here the Sturm sequences 
for the MD toggle, DD toggle, and single-operator positive feedback circuit. For simplicity of notation 
we have used x as our variable in the Sturm polynomials below, rather than the P\ eq used in the main 
text. 

The Sturm polynomials associated with the MD toggle equilibrium polynomial are: 



- 2 
lot 



fo(x) = {x-X Uot )(x + iy + xX 2 
h{x) = X 2tot + (x + l)(3x - 2X ltot + 1) 



f 2 (x) = -(2(x + l)(X ltot + l) 2 - (6x + X Uot - 2)X 2tot ) 

9X 2 2 tot ( - AX 2 l t + (X ltot (X Uot + 20) - 8)X 2 2 tot - i(X ltot + 1) ; 



-2 

9X 



4((X ltot + ir~3X 2 2 tot ) 2 



The DD Sturm polynomials are: 

fo(x) = (x 2 + lf{x - X Uot ) + xX 2tot 
h(x) = (x 2 + l)(5x 2 - 4xX ltot + 1) + X 2tot 

hix) = ^(4(x 2 + l)(x(X 1 2 tot ~5) + 6X ltot ) -X 2 tot (20x + X ltot )) 

f 3 (x) = 1 (r 2 2 tot (2X 2 tot - 4 + 4x 2 (X 2 tot - 5) - 3xX ltot (3 + Xx tot )) ~ 4(x 2 + l){X\] ot + l) 2 ) 

U(x) = -4(XxL - 5) 2 X 2 6 tot (20x + X ltot ) - {X 2 tot - bfX^Jx^X^ + 35Xi* ot - 64) - 2X ltot (XiL - 62)) 

+ lG(Xl 4 tot - AX 2 tot - b) 2 X 2tot {X ltot - x) 
f 5 (x) = - (256x1 - 3X^t t (9X 2 l t + 32X 2 L - 256) - 96^L(^L + 29X 2 1 - 8) + 256(X 2 1 + l) 3 



x 25((XiL + l) 2 + (Xllt - 5)F 2 L)' 



where 



93 = ^(^L-5) 2 



94 = 100((X lt 2 ot - 5)F 2 L + {X 2 tot + l) 2 " 
g 5 = (X 2 tot - 5) 2 (l6(5iL + l) 2 + {QXiit + 35XiL - 64)X 2 1 - 80X 



For all values of X\ tot ^ \/5, the sequence consisting of the fi(x) above may be used to determine the 

number of steady states. However, when X\ tot — > \/5, the sequence terminates prematurely (since f±{x) — > 

0) and there are problematic zeroes in the denominators of fz(x) and f§(x). We thus set X± tot = \/5 in 
the equilibrium polynomial to get 

t\l q - VM\ eq + 2P? eq - 2VM\ eq + Pi eq X 2 2 tot + F\ eq - V5 = 0 , (S13) 
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and generate a second Sturm sequence to use at X\ tot 




fo(x) 
h(x) 



x 5 - V5x 4 + 2x 3 - 2V5x 2 + xX 2tot + x - V5 
5x 4 - 4V5x 3 + 6x 2 - iVEx + X 2tot + 1 
^(24^5.t 2 - 20xX 2 2 tot - VhX 2tot + 24^5) 

--^-{40V5xX 2 6 tot - 168V5xX 2 A tot ~ 288V5xX 2 l ot + 10X 2 6 fot - 28bX 2tot + 1440X 2 ' ot ) 




h{x) 



27(256X 2tot - mX 2tot ~ lbbb2X 2tot + 55296) 



40V5(5X 2tot -21X 2tot -36) 2 



The single-operator positive feedback circuit Sturm sequence contains the following polynomials: 



, 97 (4a 3 - a 2 7 2 - 18a7 2 + 4 7 4 + 277 2 ) 

f3\x) — 9 ■ 

V ' 4(a 2 - 3 7 2 ) 2 

Recall that 7 is the rescaled degradation rate constant for the A repressor and a represents the increase in 
protein production resulting from repressor dimer binding to OR2. 

S5 Sturm sequence sign combinations 

We list in the tables below the various possible sets of Sturm sequence inequalities for the DD toggle 
(Tables SI and S2) and A repressor-PRM system (Tables S3). (The MD toggle and single-operator positive 
feedback circuit each have only one combination of inequalities that can give rise to bistability.) The 
inequality sets are written in compact form as {± ± • • • ±}. 'Allowed' refers to the logical consistency of 
the inequality set (i.e., whether all inequalities can be simultaneously satisfied) and was determined using 
Mathematica. For simplicity of notation in the tables we only show + and — in the various sequence 
positions, even though there are cases where a particular Sturm polynomial may be equal to zero without 

affecting the total number of sign changes; for example, for the DD toggle, both { — I h— } and 

{ — hO — h— } give a variation of 4. (Recall that Sturm's theorem is only concerned with consecutive 
nonzero elements and zeroes are ignored.) When zeroes were valid options for the sequence, otherwise 
strict inequalities were made nonstrict (e.g., '>' — > '>'). 
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var(Tdd,0) 


signed, 0) 


Allowed? 


var(Fdd, +oo) 


signed, +oo) 


Allowed? 


5 


{- + - + -+} 


N 


2 






4 


{- + - + —} 

{- + --+-} 
{- + + -+-} 
{- + - + +-} 


N 
Y 
Y 
Y 


1 


{+ + + + +-} 
{+ + + + —} 
{+ + + -—} 
{+ + } 


N 
N 
Y 
Y 


3 






0 


{+ + + + ++} 


N 



Table SI: Sturm sequence inequality sets for the DD toggle when X\ tot ^ The signs of the first 

two polynomials are fixed at P\ eq — > 0 and P\ eq — > +oo. Neither var{Fdd->fy = 5 nor var{Fddi +oo) = 0 
represent logically consistent sets, eliminating the need to consider any sets that could in theory satisfy 
var^Fdd-, +oo) = 2 or var^dd, 0) = 3 as candidates for bistability. 



var(Tdd, 0) 


sign (7^,0) 


Allowed? 


var{Fdd, +oo) 


signed, +oo) 


Allowed? 


4 


{- + - + -} 


N 


1 








{- + - + +} 


N 








3 


{- + + -+} 


Y 


0 


{+ + + + +} 


Y 




{- + --+} 


N 









Table S2: Sturm sequence inequality sets for the DD toggle when X\ tot = y/h. The signs of the first two and 
three polynomials are fixed at P\ eq — > 0 and P\ eq — > +oo, respectively. The set with uar(J r ^,0) = 4 is not 
logically consistent, eliminating the need to consider any sets that could in theory satisfy var^dd, +oo) = 1 
as candidates for bistability. 



var(Tp RM ,0) 




sig 




Allowed? 


var{Tp RM , 


+oo) 


sign(J"p RM ,+ 


oo) 


Allowed? 


6 


{- 


+ 


+ - + -+-} 


N 


3 








{- 


+ 


H 1 h+} 


N 






{+ H h + 


++} 


N 




{- 


+ 


H h + -+} 


Y 






{+ + + 


++} 


Y 


5 


{- 


+ 


+ - + --+} 


N 


2 




{+ + 


++} 


N 




{- 


+ 


+ -- + -+} 


Y 






{+ + 


-+} 


N 


4 






1 


{+ + 


— } 


N 



Table S3: Sturm sequence inequality sets for the A repressor-PRM system. The first four Sturm polyno- 
mials have fixed signs at P eq = 0 and P eq = +oo. Neither rar^p^^) = 6 nor var (J r p HM , +oo) = 1 
represent logically consistent sets, eliminating the need to consider any sets that could in theory satisfy 
var^p^, +oo) = 3 or var^p^, 0) = 4 as candidates for bistability. 
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S6 Analytic expressions for regions of bistability 

The regions of bistability were determined by combining the valid inequality sets and reducing to a single 
pair of inequalities in the plane of the relevant variables. Reduction was done using Mathematica. 

There is only one combination of inequalities that can give rise to bistability for the MD toggle: when 
var(T m d,0) = 3 and var ( TmA, +oo) = 0. This gives the region of bistability in X\ tot ~X2tot space as 

Xu 0 t > 8 

1/ — — 2 -~ \ — 2 1/ — — 2 — \ (S14) 

- 20X ltot + X ltot - 8 - f(X ltot ) < X 2tot < - 20X ltot + X Uot - 8 + f(X Uot ) 



with f(X ltot ) = (X ltot - 8) 3 / 2 (X lto4 )V2. 

In the case of the DD toggle, there are two different Sturm sequences we need to consider depending 
on the value of X\ tot (see Section 'Sturm polynomials' above). When X\ tot ^ \/5, the Sturm sequence 
J-dd contains six polynomials, and there are three sets of inequalities with var(J-dd, 0) = 4 and two with 
var(Tdd,+oo) = 1 that are logically consistent (Table SI). When X\ tot = y/b, the sequence Tm contains 
five polynomials, and only one set of inequalities with variTdd-, 0) = 3 and one with var{Tdd^ +oo) = 0 are 
allowed (Table S2). These inequalities may be combined to give a continuous region of bistability as the 
intersection of 

Xitot > 4 

1 /„^7-4 „_^2 nd „ /„^8 ^6 ^4 „.„^2 \l/2\ ($15) 



0 < X 2tot < — ( 9X ltot + 35X Hot -64 + 3 (9X Hot + 70X Uot + 577X Uot + 640X Uot + 1024) 

and 

256fe ot + (X 2 2 tot + l) 3 ) < SXita ( 9X 2 l t + 32X 2 L - 256^ + 96XiL f X 2tot + 29X 2tot - %] . (S16) 



As with the MD toggle, the equilibrium polynomial for the single-operator positive feedback circuit also has 
a maximum possible variation of 3, which means that the circuit exhibits bistability only in the region 

a > 9 

i (a 2 + 18a -27 -(a- 9f' 2 {a - l) 1 / 2 ) < 7 2 < ^(a 2 + 18a - 27 + (a - 9) 3 / 2 (a - l) 1 / 2 ) . (S1?) 

With the exception of the DD toggle at X\ tot = \/5 (which was treated by analyzing a second Sturm 
sequence; see 'Sturm polynomials' section above), none of the potential zeroes in the Sturm polynomial 
denominators required special treatment nor did they present any problems in determining the regions of 
bistability. 

S7 Circuit Jacobians 

The Jacobian matrices for the MD toggle, DD toggle, and single-operator positive feedback circuit are: 



J, 



md 



( -kdegi ~ k cF X 2 {t) 0 0 k baal -k cR - k cF P!(t)\ 

0 ~k deg2 - 4k kF P 2 (t) 2k kR 0 k baa2 

0 2k kF P 2 {t) -fc fefl -fc nF Xi(t) ~k nR ~k nF P 2 P 2 (t) 0 

0 0 -k nF Xi(t) -k nR - k nF P 2 P 2 (t) 0 

y -k cF X 2 (t) 0 0 0 -fc cfl -fc cF Pi(t)y 



(S18) 
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-k degl -ik iF P 1 (t) 
2k iF P 1 (t) 
0 
0 
0 
0 



2k u 
-k iR — k 0j 
0 
0 
0 



>Jf 2 (t) 



rX 2 (t) 



0 
0 

leg2-4k kF P 2 (t) 
2k kF P 2 (t) 
0 
0 



0 
0 

-k kR -k^ F X 1 (t) 
-k nF X 1 (t) 
0 



0 
0 

-k nF P 2 P 2 (t) 
-k nF P 2 P 2 {t) 
0 



?P 1 P 1 (t) 



0 
0 

-k OF p 1 p 1 (t) / 



(S19) 



(-kdeg2 - ^k kF P 2 (t) 2k kR k bas2 - k w \ 

2k kF P 2 (t) -k kR -k qF X 2 (t) -k qR -k qF P 2 P 2 (t)\ (S20) 

0 -k qF X 2 (t) -k qR -k qF P 2 P 2 {t)J 

Each of the Jacobian matrices J may be transformed to Metzler matrices Jm with a similarity transforma- 
tion: Jm = P~ 1 JP. (The P matrices for the MD toggle, DD toggle, and single-operator positive feedback 
circuit Jacobians are P mc i = diag(— 1, 1, 1, —1, 1), Pdd = diag(— 1, —1, 1, 1, —1, 1), and P p f = diag(l, 1, — 1), 
respectively.) The Jacobians are also row equivalent to the identity matrix (as confirmed with Mathemat- 
ica) and thus invertible — det(J) ^ 0 for all (positive) parameters and equilibria. 



S8 Number of steady states and stability analysis 
S8.1 Preliminaries 

Definition 1 (Positive systems) A linear system x = Ax is positive if for every nonegative initial state 
the solution x(t) is nonnegative. 

The following is a well known condition for positivity [2]: 

Theorem 1 A linear system x = Ax is positive if and only if matrix A is a Metzler matrix, i.e., its 
elements satisfy: aij > 0, V(i, j) such that i ^ j. 

Since the Jacobian matrices (shown above) are, for any choice of parameters, similar to Metzler matrices 
via linear transformations, the linearizations of systems (SI), (S2), and (S3) are positive. 

The general definition of dissipativity (see, e.g., [3]) is based on the existence of compact, forward invariant 
subsets of that absorb the system trajectories. The following definition (from [4]) is equivalent and 
easier to verify: 

Definition 2 (Dissipative systems) A system x = f(x) is dissipative if its solutions are eventually 
uniformly bounded, i.e., there exists a constant k > 0 such that: 

lim supxj(i) < k. 

t— >+oo 

Systems (SI), (S2), and (S3) are dissipative. As an example, we verify the definition for the MD toggle 
model (SI). Because the total mass of each of the DNA species X± and X2 is constant, we know that 
X±(t) < X max and -X^i) < X max V t, where X max = max{Xi tot , Xitot\- The concentration of Pi can be 
upper bounded as follows: 

Pi '(t) = -k deg iPi{t) - k cF X 2 (t)-P 1 (t) + k basl Xi{t) + k cR (X 2tot - X 2 {t)) 

< -k deg iPi{t) - k cF X 2 {t)-Pi(t) + k basl X max + k cR X max 
— ~kdeg\Plif) 4" k bas \X max + k cR X max 

< -aPi{t) + b , 
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where a = kd eg i and b = {k bas \ + k cR )X max . The right hand side of the last inequality above is a linear, 
asymptotically stable system whose solution is eventually uniformly bounded (b is a finite constant). Using 
the comparison principle [5], we conclude that Pi(t) is bounded and can find a constant k that satisfies 
the definition. 

P 2 may be similarly upper bounded. We first consider the dynamics of P/(i) = P2W + 2P2P2(i), the total 
amount of unbound P2 in the system: 

P('(t) = -k deg2 P 2 (t) - 2k nF P 2 P 2 {t)-X l {t) + 2k nR {X ltot - Xi(t)) + k bas2 X 2 (t) 

< -k deg 2P2{t) ~ 2k nF P 2 P 2 (t) ■ X l (t) + {2k nR + k bas2 )X max 

< —kdeg2P2{t) + (2& n .R + k bas2 )X max . 

The dynamics of monomeric P 2 satisfy: 

P 2 \t) = -2k kF P 2 {t) 2 - k deg2 P 2 (t) + 2k kR P 2 P 2 {t) + A; 6as2 X 2 (i) 

< -kdeg2P2{t) + 2k kR P 2 P 2 {t) + k bas2 X max 

< -k deg2 P2{t) + h R (P 2 (t) + 2P 2 P 2 {t)) + fc tes2 i mflI . 

Together, we have: 

~kdeg2 k kR \ I P 2 (t) \ / A; bas2 X maa; \ 
0 J ^(t)J \(2k nR + k bas2 )X max J 

The variables P 2 / (t) and P 2 {t) are upper bounded by a linear system with eigenvalues 

Al,2 = ^ ^-^de 9 2 ± \J kj eg2 - 4:khRkde g2 \ , 

whose real part is always negative for any value of the (positive) parameters. Being upper bounded by an 
asymptotically stable linear system, the concentrations P 2 and P/ are eventually uniformly bounded. It 
follows that P2-P2 is a ls° eventually uniformly bounded, since P2P2 < P 2 ■ 

Therefore, the ODE model of the monomer-dimer toggle system is dissipative. Note that in the absence of 
degradation (kdeg2 = 0), then we cannot reach the same conclusion (the total amount of protein will grow 
unbounded). Similar proofs can be provided for systems (S2) and (S3). 




S8.2 Stability of equilibria 

Sturm's theorem applied to the polynomial equilibrium conditions for systems (SI), (S2), and (S3) reveals 
that each system admits three positive equilibria. The stability properties of these equilibria can be 
determined by degree theory [4]. 

Definition 3 (Regular equilibrium) An equilibrium point x of system x = f(x) is regular if det(J \x)) ^ 
0 (in other words, J{x) must be invertible; alternatively, J{x) must not have eigenvalues at the origin). 

Definition 4 (Index of an equilibrium point) The index of a regular equilibrium point x is the sign 
of the determinant of —J(x): 

ind(x) = sign! det(— J(x)) ) 
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Definition 5 (Degree of a system) The degree of a dynamical system x = f(x), over a set U € R n , 
having equilibria xi, i = 1, m, is defined as: 

m 

deg(/) = ^{ind(^), s< e U, f{ Xi ) = 0} , 
i=i 

where X{ are regular equilibria. 

Theorem 2 A dissipative dynamical system x = f(x) defined on M. n has degree +1 with respect to any 
bounded open set containing all its equilibria. 

Since systems (SI), (S2), and (S3) are dissipative, by Theorem 2 all have degree +1. We further note that 
the Jacobian matrices of our systems are row equivalent to the identity matrix and thus always invertible for 
any choice of (positive) parameters and equilibria. Therefore, all equilibria are regular. To determine the 
index of each equilibrium point, we need not know the value of the equilibrium itself, since in general 

ind(x) = sign^det(— J(x))^j = sign^det(A7 — J(x))^j , with A = 0. 

Therefore, the index of an equilibrium corresponds to the sign of the constant term in the system's charac- 
teristic polynomial pj(X) = det(A/ — J(S)). For any choice of the parameters (reaction rates) in systems 
(SI), (S2), and (S3), the pj(X) have coefficients that are all positive except the constant term, which may 
be positive or negative. Thus, the sign of the constant term determines the index of the corresponding 
equilibrium. Finally, we note that the sign of the constant term in the characteristic polynomial also 
determines the stability properties of the corresponding equilibrium due to the particular structure of the 
Jacobians under consideration; we can state the following lemma [6]: 

Lemma 1 Any single equilibrium of systems (SI), (S2), and (S3) is unstable if and only if the constant 
term of the characteristic polynomial pj(X) is negative. Instability can only be driven by a simple, real 
(positive) eigenvalue. 

Proof The linearization of systems (SI), (S2), and (S3) define positive linear systems, where the Jacobians 
Jmd, Jdd, and J p f (given above as (S18), (S19), and (S20)) are similar to Metzler matrices. Therefore, 
these Jacobians always have a real dominant eigenvalue, i.e. \ max > TZe(Xi), VAj € J [7]. 

The coefficients of characteristic polynomials £>j m d(A), pjdd(A), and pj p f{X) are all real and all positive 
with the exception of the constant terms, which can be positive or negative. If the constant term of each 
Pj{\) is negative, then we know that pj(0) < 0 and it is real. In the limit A — > oo, pj(X) > 0 because all 
other coefficients are positive. Thus, there must be at least one point in the right half plane that is a root 
of pj(X). Thus, all our systems unstable, and because the various Js are similar Metzler matrices, their 
largest roots must be real. 

If the system is unstable, then the characteristic polynomial must have at least one root with positive real 
part. Ab absurdo, suppose the constant term is positive. Then instability can only occur with a pair of 
complex conjugate eigenvalues with positive real part. This is impossible because the Jacobian is a Metzler 
matrix and the dominant eigenvalue must be real. Thus, the constant term of the characteristic polynomial 
must be negative. ■ 

We can now finish our stability analysis. Our systems all have degree +1 (Theorem 2), thus when three 
equilibria are present their indices must be equal to +1, +1, and -1 so that their sum is +1 (we recall that 
all equilibria of our systems are regular). Since the index is equal to the sign of the constant term in the 
characteristic polynomial, a positive index is associated with a stable equilibrium and a negative index is 
associated with an unstable equilibrium, and we can conclude that, with three equilibria, our systems are 
bistable. Note that the unstable point does not admit local oscillatory behaviors, because local instability 
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is driven by a real eigenvalue (Lemma 1). As an alternative argument, we can also simply note that 
our systems are monotone — for any choice of parameters the Jacobians are similar to Metzler matrices, 
a property that defines a monotone system with respect to the positive orthant [8, 9] — and a monotone 
system does not admit oscillatory behaviors. 

S9 Computational support for Sturm's theorem results 

For both toggle circuits, and for each of the valid combinations of sign(J 7 , 0) and sign(J r , +oo), 1000 random 
values of X ltot and X 2to t were selected from inside and outside of the predicted bistable regions and plugged 
in to the appropriate equilibrium polynomial (Eq. (3) or (4)) which were then solved numerically. In all 
cases the number of equilibria found perfectly matched the number determined by Sturm's theorem. It is 
worth noting that the time required to test {Xi tot , X2 to t} pairs scales linearly with the number of pairs, so 
while testing a small number can be done relatively quickly, as the number of pairs becomes appreciable 
the time can be significant — up to 3 hours to test 600,000 random values of X ltot and X 2to t- 

We may also computationally check the stability of the various steady states using the circuits' Jacobian J 
and characteristic polynomial pj{\) = det (XI — J). An equilibrium point is locally stable (unstable) if and 
only all (some) of the roots of pj(X) = 0 are in the open left-half (closed right-half) of the complex plane. 
It was also recently shown that if all off-diagonal components of the Jacobian are nonnegative (that is, it 
is a Metzler matrix), or if the Jacobian may be transformed to have such a form, then any equilibrium is 
unstable if and only if the constant term of pj(X) has a sign opposite to that of all other terms in pj(X) [6]. 
We use this pj(X) constant term condition to confirm that each bistable solution set contains one solution 
representing an unstable steady state. 

The inequalities that satisfy the pj(X) constant term condition are: 

— 4 —-3 —-2 —-2 

Pleq + 4-Pl e g + 2Pi eq (X 2tot + 3) 

+ ^ e9 (4 - 2(X lfot - 2)F 2 L) ~ 2(X ltot - l)X 2tot + X 2tot + 1< 0 (S21) 

and 

(J\ eq + 2I\ eq + X 2tot + I)' - 4J\e q 0\ eq + l)X ltot X 2 2 tot < 0 (S22) 

for the MD and DD toggles, respectively. For each bistable solution set found we substituted the values 
Xitot, Xitot, an d Pi eq into Eqs. (S21) and (S22) and confirmed that only one of the three solutions satisfies 
the appropriate instability condition, supporting our claim that for the systems under consideration, if three 
roots to the equilibrium polynomial can be found, then two must be stable and one must be unstable. 
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